/***************************************************************************************
 The copyright in this software is being made available under the BSD License, included
 below. This software may be subject to other third party rights,including patent rights,
 and no such rights are granted under this license.

 Regents of the Tsinghua University, Peking University, Zhejiang University, HKUST
 Tsinghua University, Peking University, Zhejiang University, HKUST 2011

 Copyright (c) 2011, Regents of the Tsinghua University, Peking University, Zhejiang
 University, HKUST. All Rights Reserved.

 Redistribution and use in source and binary forms, with or without modification, are
 permitted provided that the following conditions are met:
   * Redistributions of source code must retain the above copyright notice,
     this list of conditions and the following disclaimer.
   * Redistributions in binary form must reproduce the above copyright notice,
     this list of conditions and the following disclaimer in the documentation
     and/or other materials provided with the distribution.
   * Neither the names of the copyright holders nor the names of its
     contributors may be used to endorse or promote products derived from this
     software without specific prior written permission.

 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
 IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
 INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
 NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
 PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
 WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
 OF SUCH DAMAGE.
***************************************************************************************/


#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include "contributors.h"
#include "global.h"
#include "ratectl.h"
#include "header.h"
#include "bbv.h"

const double THETA = 1.3636;
const int Switch = 0;

int Iprev_bits = 0;
int Pprev_bits = 0;


int Xp, Xb;
static int R;
static int Np, Nb, bits_topfield;
long T, T1;
long UpperBound1, UpperBound2, LowerBound;
double InitialDelayOffset;
const double OMEGA = 0.9;

double Wp, Wb;
int TotalPFrame;
int DuantQp;
int PDuantQp;
FILE *BitRate;
double DeltaP;

extern bbv_buffer_t *p_bbv;

/*
*************************************************************************
* Function:
* Input:
* Output:
* Return:
* Attention:
*************************************************************************
*/
void init_seq_rc()
{
    double L1, L2, L3, bpp;
    int qp;
    int i;

    Xp = 0;
    Xb = 0;

    bit_rate = input->bit_rate;
    frame_rate = ( int )( ( float )( img->frame_rate * ( input->succe_b + 1 ) ) / ( float )( input->jumpd + 1 ) );
    PreviousBit_Rate = bit_rate;



    img->total_mb_num = img->height * img->width / 256;
    if ( input->basic_unit > img->total_mb_num )
        input->basic_unit = img->total_mb_num;
    if ( input->basic_unit < img->total_mb_num )
        TotalNumberofBasicUnit = img->total_mb_num / input->basic_unit;

    MINVALUE = 4.0;


    BufferSize = bit_rate * 2.56;
    CurrentBufferFullness = 0;
    GOPTargetBufferLevel = CurrentBufferFullness;

    InitialDelayOffset = BufferSize * 0.8;


    m_windowSize = 0;
    MADm_windowSize = 0;
    img->NumberofCodedBFrame = 0;
    img->NumberofCodedPFrame = 0;
    img->NumberofGOP = 0;

    R = 0;

    if ( input->succe_b > 0 )
    {
        BETAP = 0.9;
    }
    else
    {
        GAMMAP = 0.5;
        BETAP = 0.5;
    }


    PPreHeader = 0;

    Pm_X1 = bit_rate * 1.0;
    Pm_X2 = 0.0;

    PMADPictureC1 = 1.0;
    PMADPictureC2 = 0.0;

    for ( i = 0; i < 20; i++ )
    {
        Pm_rgQp[i] = 0;
        Pm_rgRp[i] = 0.0;
        PPictureMAD[i] = 0.0;
    }
    PPictureMAD[20] = 0.0;

    PDuantQp = 2;


    PAveHeaderBits1 = 0;
    PAveHeaderBits3 = 0;
    if ( TotalNumberofBasicUnit >= 9 )
        DDquant = 1;
    else
        DDquant = 2;

    MBPerRow = img->width / 16;
    RC_MAX_QUANT = 63;
    RC_MIN_QUANT = 0;


    bpp = 1.0 * bit_rate / ( frame_rate * img->width * img->height );
    if ( img->width == 176 )
    {
        L1 = 0.1;
        L2 = 0.3;
        L3 = 0.6;
    }
    else if ( img->width == 352 )
    {
        L1 = 0.2;
        L2 = 0.6;
        L3 = 1.2;
    }
    else
    {
        L1 = 0.6;
        L2 = 1.4;
        L3 = 2.4;
    }

    if ( input->init_qp == 0 )
    {
        if ( bpp <= L1 )
            qp = 35;
        else if ( bpp <= L2 )
            qp = 25;
        else if ( bpp <= L3 )
            qp = 20;
        else
            qp = 10;
        input->init_qp = qp;
    }
}

/*
*************************************************************************
* Function:
* Input:
* Output:
* Return:
* Attention:
*************************************************************************
*/
void rc_init_gop( int np, int nb )
{
    int AllocatedBits;
    int GOPDquant;

    LowerBound = ( long )( R + bit_rate / frame_rate );
    UpperBound1 = ( long )( R + InitialDelayOffset );

    AllocatedBits = ( int )floor( ( 1 + np + nb ) * bit_rate / frame_rate + 0.5 );
    R += AllocatedBits;
    Np = np;
    Nb = nb;

    GOPOverdue = FALSE;

    img->IFLAG = 1;

    TotalPFrame = np;
    img->NumberofGOP++;
    if ( img->NumberofGOP == 1 )
    {
        MyInitialQp = input->init_qp;
        PreviousQp2 = MyInitialQp - 1;
        QPLastGOP = MyInitialQp;
    }
    else
    {
        PAverageQp = ( int )( 1.0 * img->TotalQpforPPicture / img->NumberofPPicture + 0.5 );

        GOPDquant = ( int )( 0.5 + 1.0 * ( np + nb + 1 ) / 15 );
        if ( GOPDquant > 2 )
            GOPDquant = 2;

        PAverageQp -= GOPDquant;

        if ( PAverageQp > ( QPLastPFrame - 2 ) )
            PAverageQp--;
        PAverageQp = RFM_MAX( QPLastGOP - 2, PAverageQp );
        PAverageQp = RFM_MIN( QPLastGOP + 2, PAverageQp );
        PAverageQp = RFM_MIN( RC_MAX_QUANT, PAverageQp );
        PAverageQp = RFM_MAX( RC_MIN_QUANT, PAverageQp );

        MyInitialQp = PAverageQp;
        QPLastGOP = MyInitialQp;
        Pm_Qp = PAverageQp;
        PAveFrameQP = PAverageQp;
        PreviousQp1 = PreviousQp2;
        PreviousQp2 = MyInitialQp - 1;
    }

    img->TotalQpforPPicture = 0;
    img->NumberofPPicture = 0;
    NumberofBFrames = 0;
}

/*
*************************************************************************
* Function:
* Input:
* Output:
* Return:
* Attention:
*************************************************************************
*/
void rc_init_pict( int fieldpic, int topfield, int targetcomputation )
{
    int i;


    img->coded_mb = 0;

    if ( input->channel_type == 1 )
    {
        if ( img->NumberofCodedPFrame == 58 )
            bit_rate *= 1.5;
        else if ( img->NumberofCodedPFrame == 59 )
            PreviousBit_Rate = bit_rate;
    }


    if ( ( fieldpic || topfield ) && targetcomputation )
    {
        switch ( img->type )
        {
            case P_IMG:
                if ( PreviousBit_Rate != bit_rate )
                    R += ( int )floor( ( bit_rate - PreviousBit_Rate ) * ( Np + Nb ) / frame_rate + 0.5 );

                if ( img->basic_unit == img->total_mb_num )
                {
                    if ( img->NumberofPPicture == 1 )
                    {
                        TargetBufferLevel = CurrentBufferFullness;
                        DeltaP = ( CurrentBufferFullness - GOPTargetBufferLevel ) / ( TotalPFrame - 1 );
                        TargetBufferLevel -= DeltaP;
                    }
                    else if ( img->NumberofPPicture > 1 )
                        TargetBufferLevel -= DeltaP;
                }

                else
                {
                    if ( img->NumberofCodedPFrame > 0 )
                    {
                        for ( i = 0; i < TotalNumberofBasicUnit; i++ )
                            BUPFMAD[i] = BUCFMAD[i];
                    }

                    if ( img->NumberofGOP == 1 )
                    {
                        if ( img->NumberofPPicture == 1 )
                        {
                            TargetBufferLevel = CurrentBufferFullness;
                            DeltaP = ( CurrentBufferFullness - GOPTargetBufferLevel ) / ( TotalPFrame - 1 );
                            TargetBufferLevel -= DeltaP;
                        }
                        else if ( img->NumberofPPicture > 1 )
                            TargetBufferLevel -= DeltaP;
                    }
                    else if ( img->NumberofGOP > 1 )
                    {
                        if ( img->NumberofPPicture == 0 )
                        {
                            TargetBufferLevel = CurrentBufferFullness;
                            DeltaP = ( CurrentBufferFullness - GOPTargetBufferLevel ) / TotalPFrame;
                            TargetBufferLevel -= DeltaP;
                        }
                        else if ( img->NumberofPPicture > 0 )
                            TargetBufferLevel -= DeltaP;
                    }
                }

                if ( img->NumberofCodedPFrame == 1 )
                    AWp = Wp;
                if ( ( img->NumberofCodedPFrame < 8 ) && ( img->NumberofCodedPFrame > 1 ) )
                    AWp = Wp * ( img->NumberofCodedPFrame - 1 ) / img->NumberofCodedPFrame + AWp / img->NumberofCodedPFrame;
                else if ( img->NumberofCodedPFrame > 1 )
                    AWp = Wp / 8 + 7 * AWp / 8;

                if ( input->succe_b > 0 )
                {
                    TargetBufferLevel += ( AWp * ( input->succe_b + 1 ) * bit_rate / ( frame_rate * ( AWp + AWb * input->succe_b ) ) - bit_rate / frame_rate );
                }

                break;

            case B_IMG:

                if ( PreviousBit_Rate != bit_rate )
                    R += ( int )floor( ( bit_rate - PreviousBit_Rate ) * ( Np + Nb ) / frame_rate + 0.5 );
                if ( ( img->NumberofCodedPFrame == 1 ) && ( img->NumberofCodedBFrame == 1 ) )
                {
                    AWp = Wp;
                    AWb = Wb;
                }
                else if ( img->NumberofCodedBFrame > 1 )
                {
                    if ( img->NumberofCodedBFrame < 8 )
                        AWb = Wb * ( img->NumberofCodedBFrame - 1 ) / img->NumberofCodedBFrame + AWb / img->NumberofCodedBFrame;
                    else
                        AWb = Wb / 8 + 7 * AWb / 8;
                }

                break;
        }

        if ( img->type == P_IMG )
        {

            if ( img->basic_unit == img->total_mb_num )
            {
                if ( img->NumberofCodedPFrame > 0 )
                {
                    T = ( long )floor( Wp * R / ( Np * Wp + Nb * Wb ) + 0.5 );

                    T1 = ( long )floor( bit_rate / frame_rate - GAMMAP * ( CurrentBufferFullness - TargetBufferLevel ) + 0.5 );
                    T1 = RFM_MAX( 0, T1 );
                    T = ( long )( floor( BETAP * T + ( 1.0 - BETAP ) * T1 + 0.5 ) );
                }
            }

            else
            {
                if ( ( img->NumberofGOP == 1 ) && ( img->NumberofCodedPFrame > 0 ) )
                {
                    T = ( int )floor( Wp * R / ( Np * Wp + Nb * Wb ) + 0.5 );
                    T1 = ( int )floor( bit_rate / frame_rate - GAMMAP * ( CurrentBufferFullness - TargetBufferLevel ) + 0.5 );
                    T1 = RFM_MAX( 0, T1 );
                    T = ( int )( floor( BETAP * T + ( 1.0 - BETAP ) * T1 + 0.5 ) );
                }
                else if ( img->NumberofGOP > 1 )
                {
                    T = ( long )floor( Wp * R / ( Np * Wp + Nb * Wb ) + 0.5 );
                    T1 = ( long )floor( bit_rate / frame_rate - GAMMAP * ( CurrentBufferFullness - TargetBufferLevel ) + 0.5 );
                    T1 = RFM_MAX( 0, T1 );
                    T = ( long )( floor( BETAP * T + ( 1.0 - BETAP ) * T1 + 0.5 ) );
                }
            }



            T = ( long )( ( 1.0 - 0.0 * input->succe_b ) * T );

            T = RFM_MAX( T, ( long )LowerBound );
            T = RFM_MIN( T, ( long )UpperBound2 );
        }
    }

    if ( fieldpic || topfield )
    {

        img->number_of_header_bits = 0;
        img->number_of_texture_bits = 0;


        if ( img->basic_unit < img->total_mb_num )
        {
            TotalFrameQP = 0;
            img->number_of_basic_unit_header_bits = 0;
            img->number_of_basic_unit_texture_bits = 0;
            img->total_mad_basic_unit = 0;
            NumberofBasicUnit = TotalNumberofBasicUnit;
        }
    }
    if ( p_bbv->bbv_mode )
    {
        img->bbv_delay = 0xFFFF;
    }
    else
    {
        img->bbv_delay = ( ushort_t )( p_bbv->bbv_delay * 90000.0 );
        if ( 0xFFFF == img->bbv_delay )
        {
            printf( "bbv_delay can't equal to 0xFFFF when bbv is 0xFFFF mode.\n" );
        }
    }
}

/*
*************************************************************************
* Function:
* Input:
* Output:
* Return:
* Attention:
*************************************************************************
*/
double calc_mad()
{
    int k, l;
    int s = 0;
    double MAD;

    for ( k = 0; k < 16; k++ )
        for ( l = 0; l < 16; l++ )
            s += abs( diffy[k][l] );

    MAD = s * 1.0 / 256;
    return MAD;
}

/*
*************************************************************************
* Function:
* Input:
* Output:
* Return:
* Attention:
*************************************************************************
*/
void rc_update_pict( int nbits )
{
    R -= nbits;
    CurrentBufferFullness += nbits - bit_rate / frame_rate;


    LowerBound += ( long )( bit_rate / frame_rate - nbits );
    UpperBound1 += ( long )( bit_rate / frame_rate - nbits );
    UpperBound2 = ( long )( OMEGA * UpperBound1 );

    update_bbv( p_bbv, nbits );

    return;
}

/*
*************************************************************************
* Function:
* Input:
* Output:
* Return:
* Attention:
*************************************************************************
*/
void rc_update_pict_frame( int nbits )
{

    int Avem_Qc;
    int X=0;


    if ( img->basic_unit == img->total_mb_num )
        X = ( int )floor( nbits * m_Qc + 0.5 );

    else
    {
        if ( img->type == P_IMG )
        {
            Avem_Qc = TotalFrameQP / TotalNumberofBasicUnit;
            X = ( int )floor( nbits * Avem_Qc + 0.5 );
        }
        else if ( img->type == B_IMG )
            X = ( int )floor( nbits * m_Qc + 0.5 );
    }

    switch ( img->type )
    {
        case P_IMG:
            Xp = X;
            Np--;
            Wp = Xp;
            Pm_Hp = img->number_of_header_bits;
            img->NumberofCodedPFrame++;
            img->NumberofPPicture++;
            break;
        case B_IMG:
            Xb = X;
            Nb--;
            Wb = Xb / THETA;

            img->NumberofCodedBFrame++;
            NumberofBFrames++;

            break;
    }
}


/*
*************************************************************************
* Function:
* Input:
* Output:
* Return:
* Attention:
*************************************************************************
*/
void setbitscount( int nbits )
{
    bits_topfield = nbits;
}


/*
*************************************************************************
* Function:
* Input:
* Output:
* Return:
* Attention:
*************************************************************************
*/
int update_quant_para( int topfield )
{
    double dtmp;
    int m_Bits;
    int BFrameNumber;
    int StepSize;
    int PAverageQP;
    int SumofBasicUnit;
    int i;


    if ( img->basic_unit == img->total_mb_num )
    {

        if ( ( topfield ) /*|| (img->field_control == 0)*/ )
        {
            if ( img->type == I_IMG )
            {
                m_Qc = MyInitialQp;
                return m_Qc;
            }
            else if ( img->type == B_IMG )
            {
                if ( input->succe_b == 1 )
                {
                    BFrameNumber = ( NumberofBFrames + 1 ) % input->succe_b;
                    if ( BFrameNumber == 0 )
                        BFrameNumber = input->succe_b;

                    if ( PreviousQp1 == PreviousQp2 )
                        m_Qc = PreviousQp1 + 2;
                    else
                        m_Qc = ( PreviousQp1 + PreviousQp2 ) / 2 + 1;
                    m_Qc = RFM_MIN( m_Qc, RC_MAX_QUANT );
                    m_Qc = RFM_MAX( RC_MIN_QUANT, m_Qc );
                }
                else
                {
                    BFrameNumber = ( NumberofBFrames + 1 ) % input->succe_b;
                    if ( BFrameNumber == 0 )
                        BFrameNumber = input->succe_b;

                    if ( ( PreviousQp2 - PreviousQp1 ) <= ( -2 * input->succe_b - 3 ) )
                        StepSize = -3;
                    else if ( ( PreviousQp2 - PreviousQp1 ) == ( -2 * input->succe_b - 2 ) )
                        StepSize = -2;
                    else if ( ( PreviousQp2 - PreviousQp1 ) == ( -2 * input->succe_b - 1 ) )
                        StepSize = -1;
                    else if ( ( PreviousQp2 - PreviousQp1 ) == ( -2 * input->succe_b ) )
                        StepSize = 0;
                    else if ( ( PreviousQp2 - PreviousQp1 ) == ( -2 * input->succe_b + 1 ) )
                        StepSize = 1;
                    else
                        StepSize = 2;

                    m_Qc = PreviousQp1 + StepSize;
                    m_Qc += RFM_MIN( 2 * ( BFrameNumber - 1 ), RFM_MAX( -2 * ( BFrameNumber - 1 ), ( BFrameNumber - 1 ) * ( PreviousQp2 - PreviousQp1 ) / ( input->succe_b - 1 ) ) );
                    m_Qc = RFM_MIN( m_Qc, RC_MAX_QUANT );
                    m_Qc = RFM_MAX( RC_MIN_QUANT, m_Qc );
                }
                return m_Qc;
            }
            else if ( ( img->type == P_IMG ) && ( img->NumberofPPicture == 0 ) )
            {
                m_Qc = MyInitialQp;
                img->TotalQpforPPicture += m_Qc;
                PreviousQp1 = PreviousQp2;
                PreviousQp2 = m_Qc;
                Pm_Qp = m_Qc;

                return m_Qc;
            }
            else
            {
                m_X1 = Pm_X1;
                m_X2 = Pm_X2;
                m_Hp = PPreHeader;
                m_Qp = Pm_Qp;
                DuantQp = PDuantQp;
                MADPictureC1 = PMADPictureC1;
                MADPictureC2 = PMADPictureC2;
                PreviousPictureMAD = PPictureMAD[0];


                CurrentFrameMAD = MADPictureC1 * PreviousPictureMAD + MADPictureC2;



                if ( T < 0 )
                {
                    m_Qc = m_Qp + DuantQp;
                    m_Qc = RFM_MIN( m_Qc, RC_MAX_QUANT );
                }
                else
                {
                    m_Bits = T - m_Hp;
                    m_Bits = RFM_MAX( m_Bits, ( int )( bit_rate / ( MINVALUE * frame_rate ) ) );
                    dtmp = CurrentFrameMAD * m_X1 * CurrentFrameMAD * m_X1 + 4 * m_X2 * CurrentFrameMAD * m_Bits;
                    if ( ( m_X2 == 0.0 ) || ( dtmp < 0 ) || ( ( sqrt( dtmp ) - m_X1 * CurrentFrameMAD ) <= 0.0 ) )
                        m_Qstep = ( float )( m_X1 * CurrentFrameMAD / ( double )m_Bits );
                    else
                        m_Qstep = ( float )( ( 2 * m_X2 * CurrentFrameMAD ) / ( sqrt( dtmp ) - m_X1 * CurrentFrameMAD ) );

                    m_Qc = Qstep2QP( m_Qstep );

                    m_Qc = RFM_MIN( m_Qp + DuantQp, m_Qc );
                    m_Qc = RFM_MIN( m_Qc, RC_MAX_QUANT );
                    m_Qc = RFM_MAX( m_Qp - DuantQp, m_Qc );
                    m_Qc = RFM_MAX( RC_MIN_QUANT, m_Qc );
                }
                return m_Qc;
            }
        }
    }

    else
    {

        if ( img->type == I_IMG )
        {
            m_Qc = MyInitialQp;
            return m_Qc;
        }

        else if ( ( img->type == P_IMG ) && ( img->IFLAG == 1 ) /*&& (img->field_control == 1)*/ )
        {
            m_Qc = MyInitialQp;
            return m_Qc;
        }
        else if ( img->type == B_IMG )
        {

            if ( ( topfield ) /*|| (img->field_control == 0)*/ )
            {
                if ( input->succe_b == 1 )
                {
                    BFrameNumber = ( NumberofBFrames + 1 ) % input->succe_b;
                    if ( BFrameNumber == 0 )
                        BFrameNumber = input->succe_b;
                    if ( PreviousQp1 == PreviousQp2 )
                        m_Qc = PreviousQp1 + 2;
                    else
                        m_Qc = ( PreviousQp1 + PreviousQp2 ) / 2 + 1;
                    m_Qc = RFM_MIN( m_Qc, RC_MAX_QUANT );
                    m_Qc = RFM_MAX( RC_MIN_QUANT, m_Qc );
                }
                else
                {
                    BFrameNumber = ( NumberofBFrames + 1 ) % input->succe_b;
                    if ( BFrameNumber == 0 )
                        BFrameNumber = input->succe_b;

                    if ( ( PreviousQp2 - PreviousQp1 ) <= ( -2 * input->succe_b - 3 ) )
                        StepSize = -3;
                    else if ( ( PreviousQp2 - PreviousQp1 ) == ( -2 * input->succe_b - 2 ) )
                        StepSize = -2;
                    else if ( ( PreviousQp2 - PreviousQp1 ) == ( -2 * input->succe_b - 1 ) )
                        StepSize = -1;
                    else if ( ( PreviousQp2 - PreviousQp1 ) == ( -2 * input->succe_b ) )
                        StepSize = 0;
                    else if ( ( PreviousQp2 - PreviousQp1 ) == ( -2 * input->succe_b + 1 ) )
                        StepSize = 1;
                    else
                        StepSize = 2;
                    m_Qc = PreviousQp1 + StepSize;
                    m_Qc += RFM_MIN( 2 * ( BFrameNumber - 1 ), RFM_MAX( -2 * ( BFrameNumber - 1 ), ( BFrameNumber - 1 ) * ( PreviousQp2 - PreviousQp1 ) / ( input->succe_b - 1 ) ) );
                    m_Qc = RFM_MIN( m_Qc, RC_MAX_QUANT );
                    m_Qc = RFM_MAX( RC_MIN_QUANT, m_Qc );
                }
                return m_Qc;
            }

            else
                return m_Qc;
        }
        else if ( img->type == P_IMG )
        {
            if ( ( img->NumberofGOP == 1 ) && ( img->NumberofPPicture == 0 ) )
            {
                if ( 1/*(img->field_control == 0) || ((img->field_control == 1) && (img->IFLAG == 0))*/ )
                {

                    m_Qc = MyInitialQp;
                    img->number_of_basic_unit_header_bits = 0;
                    img->number_of_basic_unit_texture_bits = 0;
                    NumberofBasicUnit--;

                    if ( ( !topfield ) && ( NumberofBasicUnit == 0 ) )
                    {
                        img->TotalQpforPPicture += m_Qc;
                        PreviousQp1 = PreviousQp2;
                        PreviousQp2 = m_Qc;
                        PAveFrameQP = m_Qc;
                        PAveHeaderBits3 = PAveHeaderBits2;
                    }
                    Pm_Qp = m_Qc;
                    TotalFrameQP += m_Qc;
                    return m_Qc;
                }
            }
            else
            {
                m_X1 = Pm_X1;
                m_X2 = Pm_X2;
                m_Hp = PPreHeader;
                m_Qp = Pm_Qp;
                DuantQp = PDuantQp;
                MADPictureC1 = PMADPictureC1;
                MADPictureC2 = PMADPictureC2;

                if ( 1/*img->field_control == 0*/ )
                    SumofBasicUnit = TotalNumberofBasicUnit;
                else
                    SumofBasicUnit = TotalNumberofBasicUnit / 2;


                if ( NumberofBasicUnit == SumofBasicUnit )
                {
                    if ( T <= 0 )
                    {
                        m_Qc = PAveFrameQP + 2;
                        if ( m_Qc > RC_MAX_QUANT )
                            m_Qc = RC_MAX_QUANT;
                        if ( topfield /*|| (img->field_control == 0)*/ )
                            GOPOverdue = TRUE;
                    }
                    else
                    {
                        m_Qc = PAveFrameQP;
                    }
                    TotalFrameQP += m_Qc;
                    NumberofBasicUnit--;
                    Pm_Qp = PAveFrameQP;
                    return m_Qc;
                }
                else
                {

                    TotalBasicUnitBits = img->number_of_basic_unit_header_bits + img->number_of_basic_unit_texture_bits;
                    T -= TotalBasicUnitBits;
                    img->number_of_basic_unit_header_bits = 0;
                    img->number_of_basic_unit_texture_bits = 0;
                    if ( T < 0 )
                    {
                        if ( GOPOverdue == TRUE )
                            m_Qc = m_Qp + 2;
                        else
                            m_Qc = m_Qp + DDquant;
                        m_Qc = RFM_MIN( m_Qc, RC_MAX_QUANT );
                        if ( input->basic_unit >= MBPerRow )
                            m_Qc = RFM_MIN( m_Qc, PAveFrameQP + 6 );
                        else
                            m_Qc = RFM_MIN( m_Qc, PAveFrameQP + 3 );

                        TotalFrameQP += m_Qc;
                        NumberofBasicUnit--;
                        if ( NumberofBasicUnit == 0 )
                        {
                            if ( ( !topfield ) /*|| (img->field_control == 0)*/ )
                            {

                                PAverageQP = ( int )( 1.0 * TotalFrameQP / TotalNumberofBasicUnit + 0.5 );
                                if ( img->NumberofPPicture == ( input->intra_period - 2 ) )
                                    QPLastPFrame = PAverageQP;

                                img->TotalQpforPPicture += PAverageQP;
                                if ( GOPOverdue == TRUE )
                                {
                                    PreviousQp1 = PreviousQp2 + 1;
                                    PreviousQp2 = PAverageQP;
                                }
                                else
                                {
                                    if ( ( img->NumberofPPicture == 0 ) && ( img->NumberofGOP > 1 ) )
                                    {
                                        PreviousQp1 = PreviousQp2;
                                        PreviousQp2 = PAverageQP;
                                    }
                                    else if ( img->NumberofPPicture > 0 )
                                    {
                                        PreviousQp1 = PreviousQp2 + 1;
                                        PreviousQp2 = PAverageQP;
                                    }
                                }
                                PAveFrameQP = PAverageQP;
                                PAveHeaderBits3 = PAveHeaderBits2;
                            }
                        }
                        if ( GOPOverdue == TRUE )
                            Pm_Qp = PAveFrameQP;
                        else
                            Pm_Qp = m_Qc;
                        return m_Qc;
                    }
                    else
                    {
                        CurrentFrameMAD = MADPictureC1 * BUPFMAD[TotalNumberofBasicUnit - NumberofBasicUnit] + MADPictureC2;
                        TotalBUMAD = 0;
                        for ( i = TotalNumberofBasicUnit - 1; i >= ( TotalNumberofBasicUnit - NumberofBasicUnit ); i-- )
                        {
                            CurrentBUMAD = MADPictureC1 * BUPFMAD[i] + MADPictureC2;
                            TotalBUMAD += CurrentBUMAD * CurrentBUMAD;
                        }

                        m_Bits = ( int )( T * CurrentFrameMAD * CurrentFrameMAD / TotalBUMAD );

                        m_Bits -= PAveHeaderBits2;

                        m_Bits = RFM_MAX( m_Bits, ( int )( bit_rate / ( MINVALUE * frame_rate * TotalNumberofBasicUnit ) ) );

                        dtmp = CurrentFrameMAD * m_X1 * CurrentFrameMAD * m_X1 + 4 * m_X2 * CurrentFrameMAD * m_Bits;
                        if ( ( m_X2 == 0.0 ) || ( dtmp < 0 ) || ( ( sqrt( dtmp ) - m_X1 * CurrentFrameMAD ) <= 0.0 ) )
                            m_Qstep = ( float )( m_X1 * CurrentFrameMAD / ( double )m_Bits );
                        else
                            m_Qstep = ( float )( ( 2 * m_X2 * CurrentFrameMAD ) / ( sqrt( dtmp ) - m_X1 * CurrentFrameMAD ) );

                        m_Qc = Qstep2QP( m_Qstep );
                        m_Qc = RFM_MIN( m_Qp + DDquant, m_Qc );

                        if ( input->basic_unit >= MBPerRow )
                            m_Qc = RFM_MIN( PAveFrameQP + 6, m_Qc );
                        else
                            m_Qc = RFM_MIN( PAveFrameQP + 3, m_Qc );

                        m_Qc = RFM_MIN( m_Qc, RC_MAX_QUANT );
                        m_Qc = RFM_MAX( m_Qp - DDquant, m_Qc );
                        if ( input->basic_unit >= MBPerRow )
                            m_Qc = RFM_MAX( PAveFrameQP - 6, m_Qc );
                        else
                            m_Qc = RFM_MAX( PAveFrameQP - 3, m_Qc );

                        m_Qc = RFM_MAX( RC_MIN_QUANT, m_Qc );
                        TotalFrameQP += m_Qc;
                        Pm_Qp = m_Qc;
                        NumberofBasicUnit--;
                        if ( ( NumberofBasicUnit == 0 ) && ( img->type == P_IMG ) )
                        {
                            if ( ( !topfield ) /*|| (img->field_control == 0)*/ )
                            {
                                PAverageQP = ( int )( 1.0 * TotalFrameQP / TotalNumberofBasicUnit + 0.5 );
                                if ( img->NumberofPPicture == ( input->intra_period - 2 ) )
                                    QPLastPFrame = PAverageQP;

                                img->TotalQpforPPicture += PAverageQP;
                                PreviousQp1 = PreviousQp2;
                                PreviousQp2 = PAverageQP;
                                PAveFrameQP = PAverageQP;
                                PAveHeaderBits3 = PAveHeaderBits2;
                            }
                        }
                        return m_Qc;
                    }
                }
            }
        }
    }
    return 0;
}


/*
*************************************************************************
* Function:
* Input:
* Output:
* Return:
* Attention:
*************************************************************************
*/
void update_rc_model()
{

    int n_windowSize;
    int i;
    double error[20], std = 0.0, threshold;
    int m_Nc;
    boolean_t MADModelFlag = FALSE;

    if ( img->type == P_IMG )
    {

        if ( img->basic_unit == img->total_mb_num )
        {
            CurrentFrameMAD = ComputeFrameMAD();
            m_Nc = img->NumberofCodedPFrame;
        }

        else
        {

            CurrentFrameMAD = img->total_mad_basic_unit / img->basic_unit;

            img->total_mad_basic_unit = 0;



            CodedBasicUnit = TotalNumberofBasicUnit - NumberofBasicUnit;
            if ( CodedBasicUnit > 0 )
            {
                PAveHeaderBits1 = ( int )( 1.0 * ( PAveHeaderBits1 * ( CodedBasicUnit - 1 ) + +img->number_of_basic_unit_header_bits ) / CodedBasicUnit + 0.5 );
                if ( PAveHeaderBits3 == 0 )
                    PAveHeaderBits2 = PAveHeaderBits1;
                else
                    PAveHeaderBits2 = ( int )( 1.0 * ( PAveHeaderBits1 * CodedBasicUnit + +PAveHeaderBits3 * NumberofBasicUnit ) / TotalNumberofBasicUnit + 0.5 );
            }
            BUCFMAD[TotalNumberofBasicUnit - 1 - NumberofBasicUnit] = CurrentFrameMAD;

            if ( NumberofBasicUnit != 0 )
                m_Nc = img->NumberofCodedPFrame * TotalNumberofBasicUnit + CodedBasicUnit;
            else
                m_Nc = ( img->NumberofCodedPFrame - 1 ) * TotalNumberofBasicUnit + CodedBasicUnit;

        }

        if ( m_Nc > 1 )
            MADModelFlag = TRUE;
        PPreHeader = img->number_of_header_bits;
        for ( i = 19; i > 0; i-- )
        {
            Pm_rgQp[i] = Pm_rgQp[i - 1];
            m_rgQp[i] = Pm_rgQp[i];
            Pm_rgRp[i] = Pm_rgRp[i - 1];
            m_rgRp[i] = Pm_rgRp[i];
        }
        Pm_rgQp[0] = QP2Qstep( m_Qc );

        if ( img->basic_unit == img->total_mb_num )
            Pm_rgRp[0] = img->number_of_texture_bits * 1.0 / CurrentFrameMAD;

        else
            Pm_rgRp[0] = img->number_of_basic_unit_texture_bits * 1.0 / CurrentFrameMAD;

        m_rgQp[0] = Pm_rgQp[0];
        m_rgRp[0] = Pm_rgRp[0];
        m_X1 = Pm_X1;
        m_X2 = Pm_X2;

        n_windowSize = ( CurrentFrameMAD > PreviousFrameMAD ) ? ( int )( PreviousFrameMAD / CurrentFrameMAD * 20 ) : ( int )( CurrentFrameMAD / PreviousFrameMAD * 20 );
        n_windowSize = RFM_MAX( n_windowSize, 1 );
        n_windowSize = RFM_MIN( n_windowSize, m_Nc );
        n_windowSize = RFM_MIN( n_windowSize, m_windowSize + 1 );
        n_windowSize = RFM_MIN( n_windowSize, 20 );
        m_windowSize = n_windowSize;
        for ( i = 0; i < 20; i++ )
        {
            m_rgRejected[i] = FALSE;
        }
        RCModelEstimator( n_windowSize );
        n_windowSize = m_windowSize;
        for ( i = 0; i < ( int )n_windowSize; i++ )
        {
            error[i] = m_X1 / m_rgQp[i] + m_X2 / ( m_rgQp[i] * m_rgQp[i] ) - m_rgRp[i];
            std += error[i] * error[i];
        }
        threshold = ( n_windowSize == 2 ) ? 0 : sqrt( std / n_windowSize );
        for ( i = 0; i < ( int )n_windowSize; i++ )
        {
            if ( fabs( error[i] ) > threshold )
                m_rgRejected[i] = TRUE;
        }
        m_rgRejected[0] = FALSE;

        RCModelEstimator( n_windowSize );

        if ( MADModelFlag )
            updateMADModel();
        else if ( img->type == P_IMG )
            PPictureMAD[0] = CurrentFrameMAD;
    }
}


/*
*************************************************************************
* Function:
* Input:
* Output:
* Return:
* Attention:
*************************************************************************
*/
void RCModelEstimator( int n_windowSize )
{
    int n_realSize = n_windowSize;
    int i;
    double oneSampleQ=0;
    double a00 = 0.0, a01 = 0.0, a10 = 0.0, a11 = 0.0, b0 = 0.0, b1 = 0.0;
    double MatrixValue;
    boolean_t estimateX2 = FALSE;

    for ( i = 0; i < n_windowSize; i++ )
    {
        if ( m_rgRejected[i] )
            n_realSize--;
    }

    m_X1 = m_X2 = 0.0;

    for ( i = 0; i < n_windowSize; i++ )
    {
        if ( !m_rgRejected[i] )
            oneSampleQ = m_rgQp[i];
    }

    for ( i = 0; i < n_windowSize; i++ )
    {
        if ( ( m_rgQp[i] != oneSampleQ ) && !m_rgRejected[i] )
            estimateX2 = TRUE;
        if ( !m_rgRejected[i] )
            m_X1 += ( m_rgQp[i] * m_rgRp[i] ) / n_realSize;
    }

    if ( ( n_realSize >= 1 ) && estimateX2 )
    {

        for ( i = 0; i < n_windowSize; i++ )
        {
            if ( !m_rgRejected[i] )
            {
                a00 = a00 + 1.0;
                a01 += 1.0 / m_rgQp[i];
                a10 = a01;
                a11 += 1.0 / ( m_rgQp[i] * m_rgQp[i] );
                b0 += m_rgQp[i] * m_rgRp[i];
                b1 += m_rgRp[i];
            }
        }
        MatrixValue = a00 * a11 - a01 * a10;
        if ( fabs( MatrixValue ) > 0.000001 )
        {
            m_X1 = ( b0 * a11 - b1 * a01 ) / MatrixValue;
            m_X2 = ( b1 * a00 - b0 * a10 ) / MatrixValue;
        }
        else
        {
            m_X1 = b0 / a00;
            m_X2 = 0.0;
        }

    }
    if ( img->type == P_IMG )
    {
        Pm_X1 = m_X1;
        Pm_X2 = m_X2;
    }
}


/*
*************************************************************************
* Function:
* Input:
* Output:
* Return:
* Attention:
*************************************************************************
*/
double ComputeFrameMAD()
{
    double TotalMAD;
    int i;
    TotalMAD = 0.0;
    for ( i = 0; i < img->total_mb_num; i++ )
        TotalMAD += img->mad_of_mb[i];
    TotalMAD /= img->total_mb_num;
    return TotalMAD;
}


/*
*************************************************************************
* Function:
* Input:
* Output:
* Return:
* Attention:
*************************************************************************
*/
void updateMADModel()
{

    int n_windowSize;
    int i;
    double error[20], std = 0.0, threshold;
    int m_Nc=0;

    if ( img->NumberofCodedPFrame > 0 )
    {

        if ( img->type == P_IMG )
        {

            if ( img->basic_unit == img->total_mb_num )
                m_Nc = img->NumberofCodedPFrame;

            else
                m_Nc = img->NumberofCodedPFrame * TotalNumberofBasicUnit + CodedBasicUnit;
            for ( i = 19; i > 0; i-- )
            {
                PPictureMAD[i] = PPictureMAD[i - 1];
                PictureMAD[i] = PPictureMAD[i];
                ReferenceMAD[i] = ReferenceMAD[i - 1];
            }
            PPictureMAD[0] = CurrentFrameMAD;
            PictureMAD[0] = PPictureMAD[0];
            if ( img->basic_unit == img->total_mb_num )
                ReferenceMAD[0] = PictureMAD[1];
            else
            {
                ReferenceMAD[0] = BUPFMAD[TotalNumberofBasicUnit - 1 - NumberofBasicUnit];
            }
            MADPictureC1 = PMADPictureC1;
            MADPictureC2 = PMADPictureC2;
        }
        n_windowSize = ( CurrentFrameMAD > PreviousFrameMAD ) ? ( int )( PreviousFrameMAD / CurrentFrameMAD * 20 ) : ( int )( CurrentFrameMAD / PreviousFrameMAD * 20 );
        n_windowSize = RFM_MIN( n_windowSize, ( m_Nc - 1 ) );
        n_windowSize = RFM_MAX( n_windowSize, 1 );
        n_windowSize = RFM_MIN( n_windowSize, MADm_windowSize + 1 );
        n_windowSize = RFM_MIN( 20, n_windowSize );

        MADm_windowSize = n_windowSize;

        for ( i = 0; i < 20; i++ )
        {
            PictureRejected[i] = FALSE;
        }
        if ( img->type == P_IMG )
            PreviousFrameMAD = CurrentFrameMAD;

        MADModelEstimator( n_windowSize );


        for ( i = 0; i < ( int )n_windowSize; i++ )
        {
            error[i] = MADPictureC1 * ReferenceMAD[i] + MADPictureC2 - PictureMAD[i];
            std += error[i] * error[i];
        }
        threshold = ( n_windowSize == 2 ) ? 0 : sqrt( std / n_windowSize );
        for ( i = 0; i < ( int )n_windowSize; i++ )
        {
            if ( fabs( error[i] ) > threshold )
                PictureRejected[i] = TRUE;
        }
        PictureRejected[0] = FALSE;

        MADModelEstimator( n_windowSize );
    }
}


/*
*************************************************************************
* Function:
* Input:
* Output:
* Return:
* Attention:
*************************************************************************
*/
void MADModelEstimator( int n_windowSize )
{
    int n_realSize = n_windowSize;
    int i;
    double oneSampleQ=0;
    double a00 = 0.0, a01 = 0.0, a10 = 0.0, a11 = 0.0, b0 = 0.0, b1 = 0.0;
    double MatrixValue;
    boolean_t estimateX2 = FALSE;

    for ( i = 0; i < n_windowSize; i++ )
    {
        if ( PictureRejected[i] )
            n_realSize--;
    }


    MADPictureC1 = MADPictureC2 = 0.0;

    for ( i = 0; i < n_windowSize; i++ )
    {
        if ( !PictureRejected[i] )
            oneSampleQ = PictureMAD[i];
    }
    for ( i = 0; i < n_windowSize; i++ )
    {
        if ( ( PictureMAD[i] != oneSampleQ ) && !PictureRejected[i] )
            estimateX2 = TRUE;
        if ( !PictureRejected[i] )
            MADPictureC1 += PictureMAD[i] / ( ReferenceMAD[i] * n_realSize );
    }

    if ( ( n_realSize >= 1 ) && estimateX2 )
    {
        for ( i = 0; i < n_windowSize; i++ )
        {
            if ( !PictureRejected[i] )
            {
                a00 = a00 + 1.0;
                a01 += ReferenceMAD[i];
                a10 = a01;
                a11 += ReferenceMAD[i] * ReferenceMAD[i];
                b0 += PictureMAD[i];
                b1 += PictureMAD[i] * ReferenceMAD[i];
            }
        }
        MatrixValue = a00 * a11 - a01 * a10;
        if ( fabs( MatrixValue ) > 0.000001 )
        {
            MADPictureC2 = ( b0 * a11 - b1 * a01 ) / MatrixValue;
            MADPictureC1 = ( b1 * a00 - b0 * a10 ) / MatrixValue;
        }
        else
        {
            MADPictureC1 = b0 / a01;
            MADPictureC2 = 0.0;
        }

    }
    if ( img->type == P_IMG )
    {
        PMADPictureC1 = MADPictureC1;
        PMADPictureC2 = MADPictureC2;
    }
}


/*
*************************************************************************
* Function:
* Input:
* Output:
* Return:
* Attention:
*************************************************************************
*/
double QP2Qstep( int QP )
{
    int i;
    double Qstep;
    static const double QP2QSTEP[8] = { 1.0, 1.0905, 1.189, 1.297, 1.414, 1.542, 1.682, 1.834 };

    Qstep = QP2QSTEP[QP % 8];
    for ( i = 0; i < ( QP / 8 ); i++ )
        Qstep *= 2;

    return Qstep;
}


/*
*************************************************************************
* Function:
* Input:
* Output:
* Return:
* Attention:
*************************************************************************
*/
int Qstep2QP( double Qstep )
{
    int q_per = 0, q_rem = 0;

    if ( Qstep < QP2Qstep( 0 ) )
        return 0;
    else if ( Qstep > QP2Qstep( 63 ) )
        return 63;

    while ( Qstep > QP2Qstep( 7 ) )
    {
        Qstep /= 2;
        q_per += 1;
    }

    if ( Qstep <= ( 1.0 + 1.0905 ) / 2 )
    {
        Qstep = 1.0;
        q_rem = 0;
    }
    else if ( Qstep <= ( 1.0905 + 1.1895 ) / 2 )
    {
        Qstep = 1.0905;
        q_rem = 1;
    }
    else if ( Qstep <= ( 1.189 + 1.297 ) / 2 )
    {
        Qstep = 1.189;
        q_rem = 2;
    }
    else if ( Qstep <= ( 1.297 + 1.414 ) / 2 )
    {
        Qstep = 1.297;
        q_rem = 3;
    }
    else if ( Qstep <= ( 1.414 + 1.542 ) / 2 )
    {
        Qstep = 1.414;
        q_rem = 4;
    }
    else if ( Qstep <= ( 1.542 + 1.682 ) / 2 )
    {
        Qstep = 1.414;
        q_rem = 5;
    }
    else if ( Qstep <= ( 1.682 + 1.834 ) / 2 )
    {
        Qstep = 1.682;
        q_rem = 6;
    }
    else
    {
        Qstep = 1.834;
        q_rem = 7;
    }

    return ( q_per * 8 + q_rem );
}
